

%% data

dChar = readtable('../../../raw_data/district_characteristics.csv');
e09 = readtable('../../../raw_data/election_results_2009.csv');
e12 = readtable('../../../raw_data/election_results_2012.csv');
cs09 = readtable('../../../raw_data/campaign_spending_2009.csv');
cs12 = readtable('../../../raw_data/campaign_spending_2012.csv');

% order observations by candidate-coalition choice
d0 = find(e12.candCM==0); N0 = size(d0,1);
d1 = find(e12.candCM==1); N1 = size(d1,1);
d2 = find(e12.candCM==2); N2 = size(d2,1); N = N0+N1+N2;
dChar = dChar([d0;d1;d2],:);
e09 = e09([d0;d1;d2],:); e12 = e12([d0;d1;d2],:);
cs09 = cs09([d0;d1;d2],:); cs12 = cs12([d0;d1;d2],:);
clear d0 d1 d2

% district indicator/aggregation matrix
A = [1:N,1:N,1:N0,N0+N1+(1:N2),1:N0+N1,1:N]';
A = 1 * (A==A');

% lagged vote shares
P = [e09.MP, ...
    e09.NA, ...
    e09.PVEM + 0.5 * e09.PM, ...
    e09.PRI + 0.5 * e09.PM, ...
    e09.PAN] ./ e09.registered;

% 2012 vote shares
S = [e12.MP ./ e12.registered; ...
    e12.NA ./ e12.registered; ...
    e12.PVEM(1:N0,1) ./ e12.registered(1:N0,1); ...
    (e12.PRI(N0+N1+(1:N2),1) + e12.PVEM(N0+N1+(1:N2),1) + e12.CM(N0+N1+(1:N2),1)) ./ e12.registered(N0+N1+(1:N2),1); ...
    e12.PRI(1:N0,1) ./ e12.registered(1:N0,1); ...
    (e12.PRI(N0+(1:N1),1) + e12.PVEM(N0+(1:N1),1) + e12.CM(N0+(1:N1),1)) ./ e12.registered(N0+(1:N1),1); ...
    e12.PAN ./ e12.registered];


%% baseline specification

% demographics
D = [dChar.femaleHH,dChar.over60,dChar.rural];
dnames = {'female','over60','rural'};
KD = size(D,2);

% menu fixed effects and demographics
X = [blkdiag(ones(N0,1),ones(N1,1),ones(N2,1)),D];
X = blkdiag(X,X,[blkdiag(ones(N0,1),ones(N2,1)),[D(1:N0,:);D(N0+N1+(1:N2),:)]],...
    [blkdiag(ones(N0,1),ones(N1,1)),D(1:N0+N1,:)],X);
NN = size(X,1);

% add log-lagged vote shares
PP = zeros(NN,3);
PP([1:N0,N+(1:N0),2*N+(1:N0),2*N+N0+N2+(1:N0),3*N+N0+(1:N0)],1) = reshape(log(P(1:N0,:)),5*N0,1);
PP([N0+(1:N1),N+N0+(1:N1),2*N+N0+N2+N0+(1:N1),3*N+N0+N0+(1:N1)],2) = reshape(log([P(N0+(1:N1),1:2),sum(P(N0+(1:N1),3:4),2),P(N0+(1:N1),5)]),4*N1,1);
PP([N0+N1+(1:N2),N+N0+N1+(1:N2),2*N+N0+(1:N2),3*N+N0+N0+N1+(1:N2)],3) = reshape(log([P(N0+N1+(1:N2),1:2),sum(P(N0+N1+(1:N2),3:4),2),P(N0+N1+(1:N2),5)]),4*N2,1);
X = [X,sum(PP,2)];

% (endogenous) campaign spending
C = [cs12.MP;...
    cs12.NA;...
    cs12.PVEM(1:N0,1);cs12.CM(N0+N1+(1:N2),1);...
    cs12.PRI(1:N0,1);cs12.CM(N0+(1:N1),1);...
    cs12.PAN];
XC = [X,C,C.^2];
K = size(XC,2);


%% baseline specification: BLP

load('../../voting_stage.mat','xi0','rBblp')

% neighbors' lagged spending, internet availability, and district travel time as instruments
nb = neighb([cs09.MP, ...
    cs09.NA, ...
    [cs09.PVEM(1:N0,1) + cs09.PM(1:N0,1); zeros(N1,1); cs09.PRI(N0+N1+(1:N2),1) + cs09.PVEM(N0+N1+(1:N2),1) + cs09.PM(N0+N1+(1:N2),1)], ...
    [cs09.PRI(1:N0,1) + cs09.PM(1:N0,1); cs09.PRI(N0+(1:N1),1) + cs09.PVEM(N0+(1:N1),1) + cs09.PM(N0+(1:N1),1); zeros(N2,1)], ...
    cs09.PAN]);
nb = [nb(:,1);nb(:,2);nb(1:N0,3);nb(N0+N1+(1:N2),3);nb(1:N0+N1,4);nb(:,5)];
itnt = dChar.internet;
itnt = [itnt;itnt;itnt(1:N0);itnt(N0+N1+(1:N2));itnt(1:N0+N1);itnt];
tvltime = dChar.travel_time;
tvltime = [tvltime;tvltime;tvltime(1:N0);tvltime(N0+N1+(1:N2));tvltime(1:N0+N1);tvltime];
Z0 = [X,itnt,tvltime,nb,nb.^2];

% Gandhi-Houde nonlinear instruments
Z0 = GH_Instr([X(:,end),C,C.^2],Z0,A);
KZ = size(Z0,2);

% SGI nodes and weights
[nu,wnu] = nwspgr('KPN',2,6);

% BLP
W = Z0'*diag(xi0.*xi0)*Z0; % GMM weighting matrix
Bblp = rBblp(:,1);
JPattern = [ones(length(S),length(Bblp)),A,sparse(length(S),KZ)];
HPattern = sparse(blkdiag([ones(length(Bblp),length(Bblp)+length(S));ones(length(S),length(Bblp)),A],ones(KZ)));
optVS = optimset('Display','off','JacobPattern',JPattern,'HessPattern',HPattern,...
    'HessFcn',@(BB,lambda)HessL(BB,lambda,XC,KD,A,[N0,N1,N2],W,nu,wnu));
[B,~,exit] = knitro_nlp(@(BB)Q_BLP(BB,length(Bblp),length(S),W),[Bblp;xi0;Z0'*xi0],[],[],...
    [sparse(KZ,length(Bblp)),Z0',-speye(KZ)],sparse(KZ,1),...
    [-Inf(length(Bblp)-2,1);zeros(2,1);-Inf(length(S)+KZ,1)],Inf(length(Bblp)+length(S)+KZ,1),...
    @(BB)Shares(BB,S,XC,KD,A,[N0,N1,N2],nu,wnu),[],optVS,'knitro.opt');
xi = B(length(Bblp)+(1:length(S)),1);

% covariance matrix
[~,~,~,Ghat] = Shares(B,S,XC,KD,A,[N0,N1,N2],nu,wnu);
Ghat = Ghat';
Ghat = - [XC,Ghat(:,K+2+(1:length(S)))\Ghat(:,K+(1:2))];
Ghat = Z0' * Ghat;
VBblp = inv(Ghat'*((Z0'*diag(xi.*xi)*Z0)\Ghat)); % robust variance
seB1 = sqrt(diag(VBblp)); % standard errors
Bblp = B(1:length(Bblp),1);
rBblp = [Bblp,seB1,2*(1-normcdf(abs(Bblp./seB1)))];

disp(' ')
disp('for Table D1 (column (V)):')
disp(print_results(round(rBblp,3),dnames,[],0,1,0))


%% region FE

% demographics
D = [dChar.region<=2,dChar.region==3,dChar.femaleHH,dChar.over60,dChar.rural];
rfe = {'north','southeast'};
KD = size(D,2);

% menu fixed effects and demographics
X = [blkdiag(ones(N0,1),ones(N1,1),ones(N2,1)),D];
X = blkdiag(X,X,[blkdiag(ones(N0,1),ones(N2,1)),[D(1:N0,:);D(N0+N1+(1:N2),:)]],...
    [blkdiag(ones(N0,1),ones(N1,1)),D(1:N0+N1,:)],X);

% add log-lagged vote shares
X = [X,sum(PP,2)];

% (endogenous) campaign spending
XC = [X,C,C.^2];
K = size(XC,2);


%% region FE: BLP

load('../../voting_stage.mat','xi0rfe','rBblp_RFE')

% neighbors' lagged spending, internet availability, and travel time as instruments
Z0rfe = [X,itnt,tvltime,nb,nb.^2];

% Gandhi-Houde nonlinear instruments
Z0rfe = GH_Instr([X(:,end),C,C.^2],Z0rfe,A);
KZ = size(Z0rfe,2);

% BLP
W = Z0rfe'*diag(xi0rfe.*xi0rfe)*Z0rfe; % GMM weighting matrix
Bblp = rBblp_RFE(:,1);
JPattern = [ones(length(S),length(Bblp)),A,sparse(length(S),KZ)];
HPattern = sparse(blkdiag([ones(length(Bblp),length(Bblp)+length(S));ones(length(S),length(Bblp)),A],ones(KZ)));
optVS = optimset('Display','off','JacobPattern',JPattern,'HessPattern',HPattern,...
    'HessFcn',@(BB,lambda)HessL(BB,lambda,XC,KD,A,[N0,N1,N2],W,nu,wnu));
[B,~,exit_RFE] = knitro_nlp(@(BB)Q_BLP(BB,length(Bblp),length(S),W),[Bblp;xi0rfe;Z0rfe'*xi0rfe],[],[],...
    [sparse(KZ,length(Bblp)),Z0rfe',-speye(KZ)],sparse(KZ,1),...
    [-Inf(length(Bblp)-2,1);zeros(2,1);-Inf(length(S)+KZ,1)],Inf(length(Bblp)+length(S)+KZ,1),...
    @(BB)Shares(BB,S,XC,KD,A,[N0,N1,N2],nu,wnu),[],optVS,'knitro.opt');
xi = B(length(Bblp)+(1:length(S)),1);

% covariance matrix
[~,~,~,Ghat] = Shares(B,S,XC,KD,A,[N0,N1,N2],nu,wnu);
Ghat = Ghat';
Ghat = - [XC,Ghat(:,K+2+(1:length(S)))\Ghat(:,K+(1:2))];
Ghat = Z0rfe' * Ghat;
VBblp = inv(Ghat'*((Z0rfe'*diag(xi.*xi)*Z0rfe)\Ghat)); % robust variance
seB1 = sqrt(diag(VBblp)); % standard errors
Bblp = B(1:length(Bblp),1);
rBblp_RFE = [Bblp,seB1,2*(1-normcdf(abs(Bblp./seB1)))];

disp(' ')
disp('for Table D1 (column (X)):')
disp(print_results(round(rBblp_RFE,3),dnames,rfe,0,1,0))


%%

save('voting_stage_IVs.mat')



